%%time
import sys, os
# get current directory
path = os.getcwd()
# get parent directory
parent_directory = os.path.sep.join(path.split(os.path.sep)[:-4])
# add utils folder to current working path
sys.path.append(parent_directory+"/subfunctions/utils")
# add integration folder to current working path
sys.path.append(parent_directory+"/subfunctions/integration")
# add HyperbolicsLCS folder to current working path
sys.path.append(parent_directory+"/demos/AdvectiveBarriers/HyperbolicLCS")
Wall time: 0 ns
In the following notebok we extract Hyperbolic LCS from tensorlines of the eigenvector field associated to the minimum/maximum eigenvalue of the Cauchy Green strain tensor from AVISO data. The notebook is structured as follows:
Hyperbolic LCS from forward computation:
Compute gradient of flow map $ \mathbf{\nabla F}_{t_0}^{t_N}(\mathbf{x}_0) $ over meshgrid using an auxiliary grid.
Compute Cauchy Green (CG) strain tensor as:
\begin{equation} \mathbf{C}_{t_0}^{t_N}(\mathbf{x}_0) = \mathbf{\nabla F}_{t_0}^{t_N}(\mathbf{x}_0)^T\mathbf{\nabla F}_{t_0}^{t_N}(\mathbf{x}_0). \end{equation}
Compute maximum eigenvalues $ \lambda_2(\mathbf{x}_0) $ of $ \mathbf{C}_{t_0}^{t_N}(\mathbf{x}_0) $.
Compute repelling LCS (extracted over the time-interval $ [t_0, t_N] $) at time $ t_0 $ as tensorlines tangent to the eigenvector field $ \mathbf{\xi_1(\mathbf{x})} $ associated to the minimum eigenvalue $ \lambda_1(\mathbf{x}) $ of $ \mathbf{C}_{t_0}^{t_N}(\mathbf{x}_0) $:
\begin{equation} \dfrac{d}{ds} \mathbf{x}(s) = \mathbf{\xi_1(\mathbf{x})} \end{equation}
For repelling LCS, the tensorlines are launched from local maxima of the $ \lambda_2(\mathbf{x}) $-field.
Compute attracting LCS (extracted over the time-interval $ [t_0, t_N] $) at time $ t_0 $ as tensorlines tangent to the eigenvector field $ \mathbf{\xi_2(\mathbf{x})} $ associated to the maximum eigenvalue $ \lambda_2(\mathbf{x}) $ of $ \mathbf{C}_{t_0}^{t_N}(\mathbf{x}_0) $:
\begin{equation} \dfrac{d}{ds} \mathbf{x}(s) = \mathbf{\xi_2(\mathbf{x})} \end{equation}
For attracting LCS, the tensorlines are launched from local minima of the $ \lambda_1(\mathbf{x}) $-field.
In case of an incompressible flow (as here) local maxima of the $ \lambda_2(\mathbf{x}) $-field coincide with local minima of the $ \lambda_1({\mathbf{x}}) $-field as $ \lambda_2(\mathbf{x}) = \dfrac{1}{\lambda_1(\mathbf{x})} $.
Note that we compute repelling and attracting LCS from a single forward computation. Similarly repelling and attracting LCS can be computed from a single backward computation. Due to the unavoidable orientational discontinuities of the eigenvector-field, the integration has to be performed manually by checking the orientation of the eigenvector at every integration-step. If the orientation is not equivalent to the orientation of the eigenvector at the previous integration step, then the eigenvector has to be flipped before integrating.
%%time
import scipy.io as sio
#Import velocity data from file in data-folder
mat_file = sio.loadmat('../../../../data/Aviso/AVISO.mat')
U = mat_file['u']
V = mat_file['v']
x = mat_file['x']
y = mat_file['y']
time_data = mat_file['t']
Wall time: 108 ms
Here we define the computational parameters and the data.
import numpy as np
# number of cores to be used for parallel computing
Ncores = 8
# time resolution of data
dt_data = time_data[0, 1]-time_data[0,0]
# periodic boundary conditions
periodic_x = False
periodic_y = False
periodic = [periodic_x, periodic_y]
# unsteady velocity field
bool_unsteady = True
# defined domain
defined_domain = np.isfinite(U[:,:,0]).astype(int)
## compute meshgrid of dataset
X, Y = np.meshgrid(x, y)
## resolution of meshgrid
dx_data = X[0,1]-X[0,0]
dy_data = Y[1,0]-Y[0,0]
delta = [dx_data, dy_data]
Here we define the spatio-temporal domain over which to consider the dynamical system.
%%time
# Initial time (in days)
t0 = 0
# Final time (in days)
tN = 25
# time step-size (in days)
dt = .1
# NOTE: In case of computing the backward FTLE field tN < t0 and dt < 0!!!!!
time = np.arange(t0, tN+dt, dt)
# length of time interval (in days)
lenT = abs(tN-t0)
# longitudinal and latitudinal boundaries (in degrees)
xmin = -3
xmax = 1
ymin = -32
ymax = -24
# make sure that domain is in the data
assert (xmax <= np.max(X) and xmin >= np.min(X) and ymin >= np.min(Y) and ymax <= np.max(Y) and t0 >= np.min(time_data) and tN <= np.max(time_data)),"The domains you are chooising are outside the domain of the data!!!!! --> redefine spatial/temporal domain"
# spacing of meshgrid (in degrees)
dx = 0.025
dy = 0.025
x_domain = np.arange(xmin, xmax + dx, dx)
y_domain = np.arange(ymin, ymax + dy, dy)
X_domain, Y_domain = np.meshgrid(x_domain, y_domain)
Wall time: 969 µs
In order to evaluate the velocity field at arbitrary locations and times, we interpolate the discrete velocity data. The interpolation with respect to time is always linear. The interpolation with respect to space can be chosen to be "cubic" or "linear".
%%time
# Import interpolation function for unsteady flow field
from ipynb.fs.defs.Interpolant import interpolant_unsteady
# set nan values to zero so that we can apply interpolant. Interpolant does not work if the array contains nan values
U[np.isnan(U)] = 0
V[np.isnan(V)] = 0
# Interpolate velocity data using cubic spatial interpolation
Interpolant = interpolant_unsteady(X, Y, U, V, time, method = "cubic")
Wall time: 85.7 ms
The Cauchy Green strain tensor $ C_{t_0}^t(\mathbf{x}) $ is computed by using an auxiliary meshgrid. 'aux_grid' specifies the ratio between the auxiliary grid and the original meshgrid. This parameter is generally chosen to be between $ [\dfrac{1}{5}, \dfrac{1}{10}] $. The computations are parallelized.
%%time
# Import package for progress bar
from tqdm.notebook import tqdm
# Import package for parallel computing
from joblib import Parallel, delayed
# Import package which checks particle location
from ipynb.fs.defs.check_location import check_location
# Import gradient of velocity function
from ipynb.fs.defs.gradient_flowmap import gradient_flowmap
# Import Rate of Strain function
from ipynb.fs.defs.CauchyGreen import CauchyGreen
# Define ratio of auxiliary grid spacing vs original grid_spacing
aux_grid_ratio = .2 # [1/5, 1/10]
aux_grid = [np.around(aux_grid_ratio*(X_domain[0, 1]-X_domain[0, 0]), 5), np.around(aux_grid_ratio*(Y_domain[1, 0]-Y_domain[0, 0]), 5)]
# parallellized computation
def parallel_C(i):
C_parallel = np.zeros((X_domain.shape[1], 2, 2))
for j in range(C_parallel.shape[0]):
x = np.array([X_domain[i,j], Y_domain[i,j]])
# only compute Cauchy Green strain tensor for trajectories starting in region where velocity field is defined
if check_location(X, Y, defined_domain, x)[0] == "IN":
# Compute gradient of flowmap
gradFmap = gradient_flowmap(time, x, X, Y, Interpolant, periodic, defined_domain, bool_unsteady, dt_data, delta, aux_grid)
# gradFmap has shape (2, 2, len(time)) --> we need gradient of flow map from t0 to tN
gradFmap_t0_tN = gradFmap[:,:,-1]
# Compute CG tensor
C_parallel[j, :, :] = CauchyGreen(gradFmap_t0_tN)
else:
C_parallel[j, :, :] = np.nan
return C_parallel
C = np.array(Parallel(n_jobs=Ncores, verbose = 0)(delayed(parallel_C)(i) for i in tqdm(range(X_domain.shape[0]))))
Wall time: 3min 41s
We now compute the eigenvalues and eigenvectors of Cauchy Green strain tensor 'C'. We additionally also compute the $ \mathrm{FTLE}_{t_0}^{t_N} $ as we seek to later compare the repelling LCS with features of the FTLE-field.
%%time
# Import eigenvalues/eigenvectors calculator
from ipynb.fs.defs.eigen import eigen
# add FTLE folder to current working path
sys.path.append(parent_directory+"/demos/AdvectiveBarriers/FTLE")
# Import FTLE calculate
from ipynb.fs.defs.FTLE import _FTLE
# maximum eigenvalues of CG-tensor
eig1 = C[:,:,0,0].copy()*np.nan
# minimum eigenvalues of CG-tensor
eig2 = C[:,:,0,0].copy()*np.nan
# eigenvectors associated to maximum eigenvalues of CG-tensor
e1 = np.zeros((C.shape[0], C.shape[1], 2))*np.nan
# eigenvectors associated to minimum eigenvalues of CG-tensor
e2 = np.zeros((C.shape[0], C.shape[1], 2))*np.nan
# FTLE-field
FTLE = C[:,:,0,0].copy()*np.nan
#iterate over meshgrid
for i in range(X_domain.shape[0]):
for j in range(Y_domain.shape[1]):
x = [X_domain[i,j], Y_domain[i, j]]
# only compute CG tensor for trajectories starting region where velocity field is defined
if check_location(X, Y, defined_domain, np.array(x))[0] == "IN":
# compute eigenvalues of rate of strain
eig1[i,j], eig2[i,j], e1[i,j,:], e2[i,j,:] = eigen(C[i,j,:,:])
# compute FTLE-field
FTLE[i, j] = _FTLE(C[i,j,:,:], lenT)
Wall time: 2.29 s
Repelling LCS can be sought among trajectories of the differential equation:
\begin{equation} \mathbf{x}'_0(s) = \xi_1(\mathbf{x}_0;t_0, t_N), \label{eq: shrinklines} \end{equation}with $ \xi_1 $ denoting the eigenvector associated to the weakest eigenvalue $ \lambda_1 $ of $ C_{t_0}^{t_N}(\mathbf{x}_0) $. The non-orientable vector field is well defiend away from tensorline singularites (points where $ \lambda_1 = \lambda_2 $). The most repelling shrinklines mark initial positions of repelling LCSs. Repelling LCSs can therefore be located as trajectories of eq. \ref{eq: shrinklines} that have locally the largest averaged $ \lambda_2(\mathbf{x}_0;t_0, t_N) $ among all neighbouring shrinklines.
%%time
from ipynb.fs.defs.tensorlines import _tensorlines
# step-size used for integration
step_size = 0.05
# threshold distance to locate local maxima in the 'eig2'
max_distance = 0.25
# maximum length of shrinkline
max_length = 5
# compute shrinklines
shrinklines = _tensorlines(X_domain, Y_domain, eig2, e1, max_distance, max_length, step_size)
Wall time: 2.51 s
Attracting LCS can be sought among trajectories of the differential equation:
\begin{equation} \mathbf{x}'_0 = \xi_2(\mathbf{x}_0;t_0, t_N), \label{eq: stretchlines} \end{equation}with $ \xi_2 $ denoting the eigenvector associated to the strongest eigenvalue $ \lambda_2 $ of $ C_{t_0}^{t_N}(\mathbf{x}_0) $. The non-orientable vector field is well defiend away from tensorline singularites (points where $ \lambda_1 = \lambda_2 $). The most attracting stretchlines mark initial positions of attracting LCSs. Attracting LCSs can therefore be located as trajectories of eq. \ref{eq: stretchlines} that have locally the smallest averaged $ \lambda_1(\mathbf{x}_0;t_0, t_N) $ among all neighbouring stretchlines.
from ipynb.fs.defs.tensorlines import _tensorlines
# step-size used for integration
step_size = 0.05
# threshold distance to locate local maxima in the 'eig2'
max_distance = 0.25
# maximum length of stretchline
max_length = 5
# compute stretchlines
stretchlines = _tensorlines(X_domain, Y_domain, eig2, e2, max_distance, max_length, step_size)
###################################### PLOT REPELLING/ATTRACTING LCS ######################################
import matplotlib.pyplot as plt
# Define figure/axes
fig = plt.figure(figsize = (16, 10), dpi = 600)
ax = plt.axes()
for i in range(len(stretchlines[0])):
ax.plot(stretchlines[0][i], stretchlines[1][i], c = 'b', linewidth = 1.5, linestyle = "dashed")
for i in range(len(shrinklines[0])):
ax.plot(shrinklines[0][i], shrinklines[1][i], c = 'r', linewidth = 1.5, linestyle = "dashed")
# add FTLE field to plot
cax = ax.contourf(X_domain, Y_domain, FTLE, levels = 400, cmap = "gist_gray")
# Set axis limits
ax.set_xlim([xmin, xmax])
ax.set_ylim([ymin, ymax])
# Set axis labels
ax.set_xlabel("long (deg)", fontsize = 16)
ax.set_ylabel("lat (deg)", fontsize = 16)
# Colorbar
cbar = fig.colorbar(cax, ticks = np.arange(0, 1, 0.1))
cbar.ax.set_ylabel(r'$ FTLE (\dfrac{1}{d}) $', rotation = 0, fontsize = 10)
# Plot legend
import matplotlib.patches as mpatches
attracting = mpatches.Patch(color='b', label='attracting LCS')
repelling = mpatches.Patch(color='r', label='repelling LCS')
plt.legend(handles=[attracting, repelling], fontsize = 16, loc = "lower left")
# Title
ax.set_title("ATTRACTING/REPELLING LCS over forward FTLE", fontsize = 20)
plt.show();
Similarly, we can obtain both repelling and attracting LCS at time $ t_N $ from a single backward computation of the backward flow map $ \mathbf{F}_{t_N}^{t_0}(\mathbf{x}_0) $. Obtaining positions of repelling LCSs at time $ t \in [t_0, t_N] $ from the forward forward computation remains problematic as it requires the forward advection of a discrete approximation of the repelling LCSs under the flow map $ \mathbf{F}_{t_0}^{t}(\mathbf{x}_0) $, generating exponentially growing errors. An analogous statement holds for attracting LCSs at time $ t_N $ by using backward advection under the backward flow map $ \mathbf{F}_{t_N}^{t}(\mathbf{x}_0) $.
Compared to the $ \mathrm{FTLE} $ field, which is only a diagnostic quantity used for visualizing possible hyperbolic LCSs, the repelling/attracting LCSs derived from the tensorline computation are precise indicators of hyperbolic LCSs. As expected repelling LCS (red solid lines) match well with ridge of the $ \mathrm{FTLE} $ field. Furthermore, due to the orthogonality of the two eigenvector fields $ \mathbf{\xi}_1 $ and $ \mathbf{\xi}_2 $, repelling LCS need always be orthogonal to attracting LCS. Intersections of repelling and attracting LCSs are referred to as generalized finite-time saddles[1, 2, 3].
[1] Farazmand, M., & Haller, G. (2012). Computing Lagrangian coherent structures from their variational theory. Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(1), 013128.
[2] Farazmand, M., & Haller, G. (2013). Attracting and repelling Lagrangian coherent structures from a single computation. Chaos: An Interdisciplinary Journal of Nonlinear Science, 23(2), 023101.
[3] Onu, K., Huhn, F., & Haller, G. (2015). LCS Tool: A computational platform for Lagrangian coherent structures. Journal of Computational Science, 7, 26-36.